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Project  2(b): 


Certified  Rapid  Solution  of  Parametrized  Partial  Differential  Equations  for 

Real-Time  Applications1 

Investigator:  Anthony  T.  Patera 

Department  of  Mechanical  Engineering 
77  Mass.  Ave.,  Room  3-266 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 

1  Introduction 

Engineering  analysis  requires  the  prediction  of  selected  “outputs”  s  relevant  to  ultimate  com¬ 
ponent  and  system  performance;  typical  outputs  include  critical  stresses  or  strains,  flowrates  or 
pressure  drops,  and  various  measures  of  temperature  and  heat  flux.  These  outputs  are  functions 
of  “inputs”  [i  that  serve  to  identify  a  particular  configuration  of  the  component  or  system;  typical 
inputs  reflect  geometry,  properties,  and  boundary  conditions  and  loads. 

In  many  cases,  the  input-output  function  is  best  articulated  as  a  (say)  linear  functional  i  of  a 
field  variable  u(fi)  that  is  the  solution  to  an  input-parametrized  partial  differential  equation  (PDE); 
typical  field  variables  and  associated  PDEs  include  temperature  and  steady/unsteady  conduction, 
displacement  and  equilibrium  elasticity/Helmholtz,  and  velocity  and  steady  incompressible  Navier- 
Stokes.  System  behavior  is  thus  described  by  an  input-output  relation  s(fj)  =  i{u(p))  the  evaluation 
of  which  requires  solution  of  the  underlying  PDE. 

Our  focus  is  on  “deployed”  systems  —  components  or  processes  in  operation  in  the  field  — 
and  associated  “Assess- Act”  scenarios.  In  the  Assess  stage  we  pursue  robust  parameter  estimation 
(inverse)  procedures  that  map  measured-observable  outputs  to  (all)  possible  system-characteristic 
and  environment-state  inputs.  In  the  subsequent  Act  stage  we  then  pursue  adaptive  design  (op¬ 
timization)  procedures  that  map  mission-objective  outputs  to  best  control- variable  inputs.  The 
computational  requirements  on  the  PDE-induced  evaluation  p  — >  s  are  formidable:  the  response 
must  be  real-time  —  we  must  “Assess- Act”  immediately ;  and  the  outputs  must  be  rigorously  cer¬ 
tified  —  we  must  “Assess- Act”  safely  and  feasibly  [23] . 

In  this  proje£o7^e“d£veiop*^j^otdi^ 

relations;  the  two  ingredients  are  reduced-basis  (RB)  approximation  [2,  8,  9,  12,  18,  22,  24,  25]  and 
a  posteriori  error  estimation  [18,  21,  27,  34,  35,  36].  We  also  consider  the  incorporation  of  these 
methods  in  the  Assess- Act  paradigm  with  particular  application  to  problem  areas  of  interest  to 
DARPA  and  GE:  in  particular,  “prognostics”  and  “accelerated  insertion  of  materials.”  In  Sections 
2-4  we  describe  the  basic  methodology  for  a  relatively  simple  class  of  problems  —  linear  elliptic 
equations;  and  in  Section  5  we  apply  this  technology  to  a  non-destructive-testing  example  —  on¬ 
line  crack  detection  —  relevant  to  ultimate  prognostics  systems.  In  Section  6  and  Section  7  we 
extend  the  methodology  to  two  classes  of  problems  relevant  to  accelerated  insertion  of  materials 
(and  materials  processing):  in  Section  6  we  consider  the  equations  of  fluid  flow;  and  in  Section  7  the 
equations  of  time-dependent  heat  transfer  (e.g.,  relevant  to  turbine  disk  quenching).  We  note  that, 
due  to  the  short  duration  of  this  project,  we  have  not  yet  applied  the  full  Assess-Act  paradigm  to 
the  models  of  Sections  6  and  7;  however,  this  will  be  pursued  in  future  work. 


1  Based  on  invited  paper  submitted  to  Volume  on  2nd  Sandia  Workshop  on  PDE- Constrained  Optimization:  Toward 
Real-Time  and  Online  PDE-Constrained  Optimization. 


2  Abstract  Statement:  Elliptic  Linear  Equations 

We  first  consider  the  “exact”  (superscript  e)  problem:  given  ^eDC  Mp,  we  evaluate  se(p)  = 
£(ue(p)),  where  ue(p)  satisfies  the  weak  form  of  our  /x-parametrized  PDE,  a(ue(p),  v\  p)  =  f(v ), 
V  v  E  Xe .  Here  p  and  V  are  the  input  and  (closed)  input  domain,  respectively;  ue(p)  is  our  field 
variable;  Xe  is  a  Hilbert  space  with  inner  product  (w,v)  and  associated  norm  \\w\\  =  yj{w^w)\  and 
a(*,  •;  p)  and  /(•),  £(•)  are  Xe-continuous  bilinear  and  linear  functionals,  respectively. 

Our  interest  here  is  in  second-order  PDEs,  and  thus  (Hq(CI))1'  C  Xe  C  (i71(0))I/;  here  C  Rd 
is  our  spatial  domain,  v  =  1  for  a  scalar  field  variable  and  v  =  d  for  a  vector  field  variable;  and 
Hl(Cl)  (respectively,  Hq(Q))  is  the  usual  Hilbert  space  of  derivative  square-integrable  functions 
(respectively,  derivative  square-integrable  functions  that  vanish  on  the  domain  boundary,  dfl)  [28]. 
The  associated  inner  product  (•,  •)  is  a  ju-independent  continuous  coercive  symmetric  bilinear  form 
over  Xe  that  perforce  induces  an  quivalent  norm  ||  •  ||. 

We  next  introduce  X  (typically,  X  C  Xe),  a  reference  finite  element  approximation  space  of 
finite  dimension  AT.  Our  reference  (or  “truth”)  finite  element  approximation  u{p)  E  X  is  then 
defined  by  a(u(p),  v;  /i)  =  f(v ),  V  v  €  X:  u(/jl)  E  X  is  a  calculable  surrogate  for  ue(p)  upon  which 
we  will  build  our  RB  approximation  and  with  respect  to  which  we  will  evaluate  the  RB  error; 
u(p)  also  serves  as  the  “classical  alternative”  relative  to  which  we  will  assess  the  efficiency  of  our 
approach.  We  assume  that  \\ue(p)  —  u(p) ||  is  suitably  small  and  hence  that  Af  is  typically  very 
large:  our  formulation  must  be  both  stable  and  efficient  as  Af  — >  oo. 

We  shall  make  two  crucial  hypotheses.  The  first  hypothesis  is  related  to  well-posedness,  and 
is  often  verified  only  a  posteriori.  We  assume  that  the  inf-sup  parameter,  (3(p)  =  iniw^x 
[a(iw,t;;/x)/(||^||||t;||)],  is  strictly  positive:  (3(p)  >  fio  >  0,  V  fi  E  V.  The  second  hypothesis  is 
related  primarily  to  numerical  efficiency,  and  is  typically  verified  a  priori.  We  assume  that  a  is 
affine  in  the  parameter  p:  a(w,  v;  p)  =  ©g(A j)a9(w,  v),  for  q  =  1, . . . ,  Q  parameter  -dependent 

functions  Qq(p)  :  V  —>  R  and  parameter- independent  continuous  bilinear  forms  aq(w,  v).  The  affine 
assumption  may  in  fact  be  relaxed  [5]. 


3  Reduced-Basis  Approximation 

The  reduced-basis  (RB)  approximation  was  first  introduced  in  the  late  1970s  in  the  context  of 
nonlinear  structural  analysis  [2,  22]  and  subsequently  abstracted  and  analyzed  [8,  25]  and  extended 
[yr  *j  2,  24]  to  a  much  larger  class  of  parametrized  partial  differential  equations.  We  first  introduce 
nested  samples  ^  i  {pi  G  Z>, . . . ,  ptv  "E^]  ■;  ‘Tagrangian” 

RB  spaces  WN  =  span  (Cn(^n)  =  u(/in)i  1  <  n  <  N},  1  <  N  <  Nmax.  Our  RB  approximation 
is  then:  given  p  E  2?,  evaluate  sx(p)  =  £(ux(p)),  where  ujsr(p)  satisfies  a{ujsf{p)^v\p)  =  f(v), 
V  v  E  Wn.  We  consider  here  only  Galerkin  projection. 

In  essence,  Wn  comprises  “snapshots”  on  the  parametrically  induced  manifold  M  =  {u(p)  \  p  E 
V }  C  X.  It  is  clear  that  M.  is  very  low- dimensional]  furthermore,  it  can  be  shown  under  our 
hypotheses  —  we  consider  the  equations  for  the  sensitivity  derivatives  and  invoke  stability  and 
continuity  —  that  M  is  very  smooth.  We  thus  anticipate  that  ux(p)  — >  u(p)  very  rapidly,  and 
hence  that  —  at  least  for  modest  P  —  we  may  choose  TV  <  AT.  Many  numerical  examples  justify  this 
expectation  (see  Sections  5,  6,  and  7);  and,  in  certain  simple  cases,  exponential  convergence  can  be 
proven  [19].  We  emphasize  that  the  deployed  context  requires  global  reduced-basis  approximations 
that  are  uniformly  (rapidly)  convergent  over  the  entire  parameter  domain  V\  proper  choice  of  the 
parameter  samples  Sn  is  thus  crucial  (see  Section  4). 


We  now  represent  ux(p)  as  uw(p)  =  where  N  =  {l,...,iV},  and  Nmax  = 

{1, . . . ,  Nmax}*  Our  RB  output  may  then  be  expressed  as  sn(p)  =  uNj(l*)&{Cj)i  where  —  we 
now  invoke  our  affine  assumption  —  the  uxj{p ),  1  <  j  <  N,  satisfy  the  N  x  N  linear  algebraic 
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system 


(i) 


E  {  E  ©%)<*%•, C i)}«JVi(M)  =  /(C<),  Vi  g  N  , 

j€N  g€Q 

where  Q  =  {1,...,Q}.  (In  practice  we  replace  the  (j,  1  <  j  <  AT,  with  a  (•, -)-orthonormalized 
system;  the  algebraic  stiffness  matrix  is  then  well-conditioned.)  It  is  clear  from  (1)  that  we  may 
pursue  an  offline-online  computational  strategy  [3,  12,  18,  27]  ideally  suited  to  the  deployed  real¬ 
time  context. 

In  the  offline  stage  —  performed  once  —  we  first  solve  for  the  Q,  V  i  G  Nmax;  we  then  form 
and  store  £(Q),  Vi  <E  Nrnax,  and  V  (i,j)  G  N^,  V  q  G  Q.  In  the  online  stage  — 

performed  many  times,  for  each  new  p  “in  the  field”  —  we  first  assemble  and  subsequently  invert 
the  (full)  N  x  N  “stiffness”  matrix  E9gq©9 (Cj > Ci)  to  obtain  the  Uffj(p),  1  <  j  <  N  —  at 
cost  0(QN2)  +  0(7V3);  we  then  evaluate  the  sum  EjsN UN to  obtain  sn(p)  —  at  cost 
O(N).  The  online  complexity  is  independent  of  J\f,  and  hence  —  given  that  N  <C  M  —  we  shall 
realize  extremely  rapid  “deployed”  response. 


4  A  Posteriori  Error  Estimation 

We  first  “presume”  P(p),  a  (to-be-constructed)  positive  lower  bound  for  the  inf-sup  parameter, 
••  j3(p)\  fflp)  >  f}(p)  >  p0  >  0,  V  p  G  V.  We  next  introduce  the  dual  norm  of  the  residual: 
eN{fj)  =  supv€JC[-R(v; /x)/||t;|| ],  where  R{v\p)  =  f(v)  -  flW/i),^),VvG  X . 

We  may  now  define  our  error  estimator,  Ax(aO  =  ^(aO/^O-Oj  and  associated  effectivity, 
Vn{^)  =  [Ajv(/z)/||it(/z)  -  u*Qi)||].  We  can  then  readily  demonstrate  [27,  36]  that 

1  <  mil*)  <  7 V/iGP,  V  TV  G  Nmax,  (2) 


where  7 (/i)  =  supwGXsupveX[a(u;,t;;/i)/(||t(;||||v||)]  is  our  continuity  “constant.”  The  left  inequality 
states  that  Aw(aO  is  a  rigorous  upper  bound  for  || u(ijl)  -  ux(/i) ||;  the  right  inequality  states  that 
A pj(fi)  is  a  (reasonably)  sharp  upper  bound. 

We  may  also  develop  bounds  for  the  error  in  the  output;  we  consider  here  the  special  “compli¬ 
ance”  case  in  which  £  =  /  and  a  is  symmetric  —  more  general  functionals  £  and  nonsymmetric  a 
require  adjoint  techniques  [27].  We  first  define  our  output  error  estimator,  A sN(fji)  = 
which  scales  as  the  square  of  the  dual  norm  of  the  residual,  ew(/i).  We  can  then  demonstrate 
[21,  27;  .36]  that  1  r  V,  V  ax\r  a  rigorous  upper 

bound  for  |s(/z)  —  sw(aOI>  we  may  further  prove  [27]  in  the  coercive  case  that  AsN(fj)/\s(fj,)  —  snO^)!  ^ 
7 —  Ax(/i)  is  a  (reasonably)  sharp  upper  bound. 

It  remains  to  develop  appropriate  constructions  and  associated  offline-online  computational 
procedures  for  the  efficient  calculation  of  £n(v>)  and  Pif1)-  To  begin,  we  consider  the  former  [18, 
21,  27]:  we  invoke  duality,  our  reduced-basis  expansion,  the  affine  parametric  dependence  of  a,  and 
linear  superposition  to  express 

e2M  =  (c,c)+  E  E0,(/‘WnM{2(c,4)+  E  E  ©v(m)Wm)(£U*;)}> 

geQneN  qr'eQn'CN 


where  C  G  X  and  Cn  G  X,  V  n  G  N,  V  q  G  <Q>  satisfy  the  parameter-independent  Poisson(-like) 
problems  (C,  v)  =  f(v),  V  v  G  X  and  (Cn,  v)  =  — a9(£„,  v),  V  v  G  X. 

An  efficient  offline-online  decomposition  may  now  be  identified.  In  the  offline  stage  —  performed 
only  pnce  —  we  first  solve  for  C  and  Cn,  Vn  G  Nmax>  V <j  G  Q:  we  then  form  and  store  the  associated 
parameter-independent  inner  products  (C,C),(C,Cn),(Cn,Cqn,),  V  (n,  n')  G  N^,  V  (q,  q')  G  Q2.  In 
the  online  stage  —  performed  many  times,  for  each  new  value  of  p  “in  the  field”  —  we  simply 
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evaluate  the  e2N(p)  sum  in  terms  of  Qq(n),  ujsrn(ii),  and  the  precomputed  inner  products  —  at  cost 
0(Q2N2).  The  online  cost  is  independent  of  N  and  —  for  Q  not  too  large  —  commensurate  with 
the  online  cost  to  evaluate  sat(/j). 

Finally,  we  turn  to  the  development  of  our  lower  bound  /?(/i)  for  the  inf-sup  “constant”  /3(/i). 
For  simplicity,  we  consider  here  the  particular  case  P  =  1,  Q  =  2,  01(/x)  =  1,  Q2(fx)  =  /i: 
a(w,v\p)  =  al(w,v)  4-  jia2(w,v)\  we  further  suppose  that  V  is  convex.  The  more  difficult  general 
case  is  considered  in  [21]  and  illustrated  in  subsequent  sections.  As  our  point  of  departure,  we  note 
that  =  infvex  v^&(v,v;m)/IMI2\  where  b(w,v\fi)  =  ( T^w.T^v ),  Vid,veX,  and  w  G  X  -> 
T^w  £  X  is  defined  as  u)  =  a(u;,  v\  /x),  V  i>  G  X. 

Next,  given  any  peD  and  constant  €  G  ]0, 1[ ,  we  introduce  t(w,  v;n\p)  =  b(w,v\JX)  +  (fi- 
Ji)[a2 (w^T^v)  +  a2fy,T^w)]  and  TP  =  {fi  G  £>  |  £(u,  v;  fi\ ~p)  >  0}.  We  may  then  define  r(pi\JL)  = 
infy€x  M) /Jh^ll2.  V/i  GF.  Our  function  r(>;/J)  enjoys  three  properties:  (i)  /?(^)  > 

r(fi\JL)  >  0,  V/i  G  P*1;  (m)  r(/i;/l)  is  concave  in  p  over  the  convex  domain  V^\  and  (in)  r(fi;Ji) 
is  “tangent”  to  (3(p)  at  /z  =  p.  (To  make  property  (Hi)  rigorous  we  must  in  general  consider 
non-smooth  analysis  and  also  possibly  a  continuous  spectrum  as  M  — >  oo.) 

We  can  now  develop  our  inf-sup  lower  bound  f3  :  V  —>  R.  We  first  require  a  sample  Ej  =  {fii  G 
£>, . . .  ,/Jj  G  V}  and  associated  set  of  polytopes  Cj  =  {'Pi  C  P^1, . . .  ,Pj  C  P^J}  that  satisfy  (a)  a 
“Positivity  Condition,”  r(//;pj)  >  e  V/iG  Pj,  1  <  j  <  J,  and  (b)  a  “Coverage  Condition,” 

V  C  we  may  then  define 


hv)  = 


max  e 


/%)* 


(3) 


(We  can  also  develop  piecewise  linear  approximations,  though  —  as  discussed  further  in  Section  5 
■ —  our  inf-sup  lower  bound  need  not  be  highly  accurate.)  It  is  readily  demonstrated  that  P(p)  has 
the  requisite  theoretical  and  computational  attributes:  (3(p)  >  j3(fi)  >  efto  >  0,  V/z  G  V;  the  online 
complexity  pt  — ►  j3(fi)  depends  only  (and  at  most  linearly)  on  J,  which  in  turn  depends  only  on  P 
and  Q  and  not  on  AT.  Properties  (i),  (ii),  and  (Hi)  permit  us  to  parlay  relatively  few  expensive 
(offline)  evaluations  into  a  very  inexpensive  global  (online)  lower  bound. 

As  an  illustrative  example  we  consider  the  Helmholtz-elasticity  crack  example  of  the  next  section 
for  n  —  (u2  G  [2.5, 5.0],  0  =  1.0,  L  =  0.2)  —  the  crack  location,  z,  and  crack  length,  L,  are  fixed,  and 
only  the  frequency  squared.  u2< ,  is — r  and  pa^  ”  0.1. 

We  find  that  a  sample  Ej= 3  suffices  to  satisfy  our  Positivity  and  Coverage  Conditions  for  c  =  0.4. 
We  present  in  Figure  1(a)  (3(fj)\  r(/z;/fy)  for  /i  G  ,  1  <  j  <  J;  and  our  lower  bound  (3).  We 
note  that  (3(fj)  is  not  concave  (or  convex)  or  even  quasi-concave,  and  hence  r(/z;  Ji)  is  a  necessary 
intermediary  in  the  construction  of  our  lower  bound. 

In  conclusion,  we  can  calculate  a  rigorous  and  sharp  upper  bourfd  for  | s(fj)  —  sat(/z)|,  AsN(p)  = 
e2N(fi)/P(fj),  with  online  complexity  independent  of  Af.  These  inexpensive  error  bounds  serve  most 
crucially  in  the  deployed  stage  - —  to  choose  optimal  JV,  to  confirm  the  desired  accuracy,  to  establish 
strict  feasibility,  and  to  control  sub-optimality.  However,  the  bounds  may  also  be  gainfully  enlisted 
in  the  pre-deployed  stage  —  to  construct  optimal  samples  Sn  [21,  36]:  Given  S°pt  =  fi*  [DO 
N  =  2, ATmax;  n*N  =  argmaxM€EF  A^_x(/x);  Sff  =  U  n*N;  END  ];  our  input  sample  EF  can 
be  large  since  the  marginal  cost  to  evaluate  AsN(fj)  is  small.  (In  contrast  to  POD  economization 
procedures  [31]  we  never  form  the  rejected  snapshots :  our  inexpensive  bound  AsN(p)  serves  as  a 
(good)  surrogate  for  the  actual  error.) 


6 


(a) 


Figure  1:  Helmholtz-elasticity  example:  (a)  Plol 
(b)  Crack  parameter  uncertainty  region  TL. 


(b) 


z 

of  for  ji  £  1  <  j  <  J;  and 


5  Assess- Act  Example:  Helmholtz-Elasticity 

We  apply  the  RB  method  here  to  a  Helmholtz-elasticity  equation  often  encountered  in  solid 
mechanics:  inverse  analyses  based  on  the  Helmholtz-elasticity  PDE  can  gainfully  serve  in  non¬ 
destructive  evaluation  (NDE)  procedures  for  crack  characterization  [10,  15,  17]  and  damage  assess¬ 
ment  [13,  16].  The  RB  method  significantly  improves  the  efficiency  of  these  inverse  procedures  — 
accelerating  the  many  evaluations  [14]  of  the  PDE  outputs. 

We  consider  a  two-dimensional  thin  plate  with  a  horizontal  crack  at  the  (say)  interface  of  two 
lamina:  the  (original)  domain  ft°(z,  L)  C  R2  is  defined  as  [0,2]  x  [0,1]  \Pq,  where  P£  =  {x\  £ 
[z  -  L/2,z  +  L/ 2],X2  =  1/2}  defines  the  idealized  crack.  The  left  surface  of  the  plate  is  secured; 
the  top  and  bottom  boundaries  are  stress-free;  and  the  right  boundary  is  subject  to  a  vertical 
oscillatory  uniform  force  of  frequency  u>.  We  model  the  plate  as  plane-stress  linear  isotropic  elastic 
with  (scaled)  density  unity,  Young’s  modulus  unity,  and  Poisson  ratio  0.25;  the  latter  determine  the 
(parameter-independent)  constitutive  tensor  Eijki ■  Our  input  is  fx  =  =  (w2,  z,  L); 

our  output  is  the  (oscillatory)  amplitude  of  the  average  vertical  displacement  on  the  right  edge  of 
the  plate. 


Thx:.  >.f''*r-“xhr,..dj£pIciC6iiiout  yP (x° *  t;\  G  ic  tlicrsfcrc  ji) ■  ■ 

f°(v),  Vu  £  where  is  a  quadratic  finite  element  truth  approximation  subspace  (of  di¬ 

mension  M  —  14,662)  ofXe(/z)  =  {n  £  (ff1(fl°(z,L)))2  |  v|x°=o  =  0};  a°(w,v,fi)  =  $Q,°(ZtL)wijEijkt 
Vk,e  ~  u2WiVi  (vij  denotes  dvi/dxj  and  repeated  physical  indices  imply  summation);  and  f°(v)  = 
f  0=2V2-  The  crack  surface  is  hence  modeled  extremely  simplistically  —  as  a  stress-free  boundary. 
(No  crack-tip  element  is  needed  as  the  output  of  interest  is  far  from  the  crack.)  The  output  s°(/x) 
is  given  by  s°(/i)  =  where  £°(v)  =  f°(v);  we  are  thus  “in  compliance.”  (For  the  damped 

example  of  Section  4  we  suitably  complexify  our  field  variable  and  space  and  replace  Eijke  with 
a  very  simple  “hysteretic”  Kelvin  model  [4]  Etjke{l  +  here  dm  is  a  material  damping 

constant.) 

We  now  map  ft°(z,  L )  via  a  continuous  piecewise-affine  transformation  to  a  fixed  domain  ft. 
This  new  problem  can  now  be  cast  precisely  in  the  desired  abstract  form  of  Section  2,  in  which 
ft,  X,  and  (w,  v)  are  independent  of  the  parameter  /i:  as  required,  all  parameter  dependence  now 
enters  through  the  bilinear  and  linear  forms.  Furthermore,  it  is  readily  demonstrated  that  our 
affine  assumption  applies  for  Q  =  10;  the  ©9(/r)  are  of  the  form  M(i)M(2)^(3)  f°r  exponents  y\  =  0 
or  1,  1/2  =  -1,0,  or  1,  and  1/3  =  -1,0,  or  1.  See  [21]  for  detail  of  the  Qq{(i),  aq(w,v),  1  <  q  <  Q, 
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and  the  “bound  conditioner”  (-,*)♦ 

We  shall  consider  two  different  models.  In  Model  I,  relevant  to  the  Assess  stage,  we  consider 
the  parameter  domain  V1  =  [3.2, 4.8]  x  VZ,L ,  where  VZ,L  =  [0.9, 1.1]  x  [0.15,0.25].  Note  that  V1 
does  not  contain  any  resonances,  and  hence  P{(i)  is  bounded  away  from  zero;  however,  u2  =  3.2 
and  cj2  =  4.8  are  quite  close  to  corresponding  natural  frequencies,  and  hence  the  problem  is 
distinctly  non-coercive.  In  Model  II,  relevant  to  the  Act  stage,  we  consider  the  parameter  domain 
pH  -  [w2  =  o]  x  £>*’L,  where  —  as  in  Model  I  —  VZ'L  =  [0.9, 1.1]  x  [0.15, 0.25].  Note  that  Model  II 
is  essentially  steady  linear  elasticity  and  thus  the  problem  is  coercive  and  relatively  easy;  we  shall 
hence  focus  our  attention  on  Model  I. 

We  first  present  basic  numerical  results.  For  our  reduced-basis  spaces  we  pursue  the  optimal 
sampling  strategy  described  in  Section  4  for  N ^  =  32  (Model  I)  and  =  6  (Model  II);  for  our 
inf-sup  lower  bound  samples  we  choose  e  =  1/5  which  yields  J1  =  84  and  J11  =  1.  We  present  in 
Table  1  A^max^ei,  VN, ave,  A^max,  and  7^ave  as  a  function  of  N  =  Nl  for  Model  I.  Here  A^max,rei  is 
the  maximum  over  Sxest  of  A^(^)/||w^max(^)||maxJ  Vn, ave  is  the  average  over  Srest  of  A;v(/i)/|| u(/i)— 
uN(fJ,)\l  A^max  rel  is  the  maximum  over  STest  of  AsN(fi)/\sNjn!lx(fi)\maxl  and  7 fN^ve  is  the  average 
over  Hxest  of  A^(/x)/|s(/x)-S7v(/i)|.  HereSTest  G  (P1)343  is  a  random  parameter  sample  of  size  343, 
II^Nmax (aO Umax  =  nia x^^HXest  ||7i;vmax  (aO  II  =  2.0775,  and  |sjvmax(/i)|max  =  max^^^Test  |sjvmax  (aO I 
0.089966.  We  observe  that  the  RB  approximation  converges  very  rapidly,  and  that  our  rigorous 
error  bounds  are  in  fact  quite  sharp.  The  effectivities  are  not  quite  0(1)  primarily  due  to  the 
relatively  crude  inf-sup  lower  bound.  (Thanks  to  the  rapid  convergence  of  RB  approximations, 
0(10)  effectivities  do  not  significantly  (adversely)  affect  efficiency.) 


N 

A,/v,max,rel 

VN}a.ve 

As 

N%  max,rel 

Vn,  ave 

10 

6.19B-01 

13.11 

8.40E-01 

22.50 

15 

5.76E-02 

13.44 

4.74E-03 

17.22 

20 

1.58E-02 

13.22 

4.50E-04 

15.44 

25 

5.69E-03 

12.57 

4.47E-05 

14.50 

30 

1.32E-03 

12.47 

2.95B-06 

14.27 

Table  1:  Numerical  results  for  Model  I. 

Turning  now  to  computational  effort  (again  for  Model  I),  for  N  =  TV1  —  25  and  any  given  /z 
(say,  4.0, 1.0, 0.2)  —  for  which  the  error  in  the  reduced-basis  output  s;v(/i)  relative  to  the  truth 
(approximation)  s(/z)  is  certifiably  less  than  A^(/z)  (say,  2.38 x  10~7)  —  the  Online  Time  to  compute 
both  sat(m)  and  Aj^(/z)  is  less  than  1/330  times  the  Time  to  directly  calculate  s(/z)  =  Z{u(fi)). 
Clearly,  the  savings  will  be  even  larger  for  problems  with  more  complex  geometry  and  solution 
structure  in  particular  in  three  space  dimensions.  Nevertheless,  even  for  our  current  very  modest 
example,  the  computational  economies  are  very  significant. 

We  now  consider  an  Assess- Act  scenario  that  illustrates  the  new  capabilities  enabled  by  rapid 
certified  input-output  evaluation  [1].  We  first  consider  the  Assess  stage  (based  on  Model  I):  given 
experimental  measurements  in  the  form  of  intervals  [s(a;2,  z*,L*)(  1  —  eexp),  $(w2,  z*,  L*)(  1  +  eexp)], 
1  <  k  <  K,  we  wish  to  determine  a  region  U  G  VZ'L  in  which  the  true  but  unknown  crack 
parameters,  (z*,  L*),  must  reside.  We  first  introduce  s^(/z)  =  ±  A^(/z),  and  recall  that  — 

thanks  to  our  bound  theorem  (2)  —  s(fi)  G  [s^(/z),  Sjv(aO]*  We  may  then  define 

K  =  {(z,L)  e  VZ'L  I  [s^(wfc,  2,  L), s^(wfc,  2,  L)\  n  [s(uil,  2*,  L*) 

(1  -  eexp),  s( cjI  z*,  L*)(l  +  eexp)]  ±  0, 1  <  k  <  K}- 
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clearly,  we  have  accommodated  both  numerical  and  experimental  error  and  uncertainty  (within 
our  model  assumptions),  and  hence  ( z*,L *)  G  71. 

In  Figure  1(b)  we  present  72.  for  if  =  2  and  (ujf  =  3.2,  u)$  =  4.8)  for  eexp  =  0.5%,  1%,  5%. 
(In  actual  practice,  we  first  find  one  point  in  7l\  we  then  conduct  a  binary  chop  at  different 
angles  to  map  out  the  boundary  of  71.)  As  expected,  as  eexp  decreases,  TZ  shrinks  towards  the 
exact  (synthetic)  value,  z*  =  1.05,  L*  =  0.17.  More  importantly,  for  any  finite  eexp,  71  rigorously 
captures  the  uncertainty  in  our  assessment  of  the  crack  parameters  without  a  priori  regularization 
hypotheses  [7].  The  crucial  new  ingredient  is  reliable  fast  evaluations  that  permit  us  to  conduct 
a  much  more  extensive  search  over  parameter  space;  for  a  given  eexp,  71  may  be  generated  online 
in  less  than  51  seconds  on  a  Pentium  1.6  GHz  laptop.  Our  search  over  possible  crack  parameters 
will  certainly  never  be  truly  exhaustive,  and  hence  there  may  be  small  undiscovered  “pockets  of 
possibility”  in  VZ'L\  however,  we  have  clearly  reduced  the  uncertainty  relative  to  more  conventional 
approaches.  (Of  course,  our  procedure  can  also  only  characterize  cracks  within  the  specified  low¬ 
dimensional  parametrization;  however,  more  general  null  hypotheses  can  be  constructed  to  detect 
model  deviation.) 

Finally,  we  consider  the  Act  stage  (based  on  Model  II).  We  presume  here  that  the  component 
must  withstand  an  in-service  steady  force  (normalized  to  unity)  such  that  the  deflection  s(0,  z*,  L*) 
in  the  “next  mission”  does  not  exceed  a  specified  value  smax.  Of  course,  in  practice,  we  will  not 
be  privy  to  ( z*,L *).  To  address  this  difficulty  we  first  define  s£  =  max(2|L)6K,  s^(0,  z,  L),  where 
s^(0,  z,L)  -  sn(0,z,L)  +  ASN(0,  z,L)\  our  corresponding  “go/no-go”  criterion  is  then  given  by 
s+  <  smax.  It  is  readily  observed  that  rigorously  accommodates  both  experimental  (crack) 
and  numerical  uncertainty  —  .s(0,  z*,L*)<s+  —  and  that  the  associated  go/no-go  discriminator 
is  hence  fail-safe.  Furthermore,  as  eexp  tends  to  zero  and  iV1  and  N 11  increase,  will  tend  to 
$(0,  z*,L*);  indeed,  for  eeXp  =  1%  and  Nl  =  25,  Nu  =  6,  [s^— s(0,  L*)]/|s(0,  z*,L*)\  =  4.73E-05. 

In  summary,  in  real-time,  we  can  both  Assess  the  current  state  of  the  crack  and  subsequently  Act 
to  ensure  the  safety  (or  optimality)  of  the  next  “sortie.” 

6  Incompressible  Navier-Stokes  Equations 

To  illustrate  the  difficulties  that  arise  in  the  treatment  of  nonlinear  problems  we  consider  a 
particular  example  [34]:  the  steady  incompressible  Navier-Stokes  equations  —  Pr(andtl)  =  0  natural 
convection  in  an  enclosure  (29,  32].  -  <' 

-  Our  ibrniuiation1  of  Section^  is  •  star  applicable '  (except  r’a  is  no  longer  =  Gr 

Grashof  number;  V  =  [1.0, 1.0E5];  ue(fi)  =  (lif  (/x) ,  t/2 (m) )  is  the  velocity  field;  fi  =  [0,4]  x  [0,1]; 
Xe  -  |  V-v  =  0};  (w,?;)  =  fQ u)i,jVitj]  a{w,v)  =  ao{w,v)  +  ±a i(w,w,v),  where  a0(w,z;)  = 

In  wiJvi,j  and  ai(w>  z>v)  =  -  Jn(wizj+wjzi)vhj  are  the  viscous  and  convective  terms,  respectively; 
f(v;p)  =  nfo(v)  =  p  fn  (l  “  \xi)  v2  is  the  buoyancy  contribution;  and  £(v)  =  2  JTq  vi(x)dx2  (for 
r0  =  {x\  =  2,X2  €  [0.5;  1]})  measures  the  flowrate.  (Note  that  the  pressure  does  not  appear 
explicitly  since  we  pose  the  problem  over  divergence-free  velocity  fields.) 

We  next  introduce  X,  a  reference  finite  element  approximation  space.  Our  reference  (or  “truth”) 
finite  element  approximation  u(p)  G  X  is  then  defined  by  a(u(p),v)  =  f(v\p),  \/v  £  X.  As  before, 
u{p)  is  a  surrogate  for  ue(p)  upon  which  we  build  our  RB  approximation,  and  relative  to  which  we 
measure  our  RB  error.  Here  X  is  the  space  (of  dimension  M  =  2,786)  of  discretely  divergence-free 
functions  associated  with  a  classical  Taylor-Hood  P2  —  Pi  finite  element  approximation  [9].  (For 
future  reference,  we  also  define  A,  the  full  Taylor-Hood  velocity  space.) 

The  derivative  of  a  plays  a  central  role:  here  da(w,  v;z)  =  ao(w,v)  4-  ai(w,z,  v)  satisfies 
a(z  +  w,v)  =  a(z,v)  +  da(w,v;z)  +  \a\(w,w,v).  It  is  readily  shown  [34]  that  da(w,v\z)  < 
j{z)\\w\\\\v\\  for  7 (z)  =  1  +  p2\\z\\\  here  p  =  v^supvGjf  IMIz,4(Q)/IMI  is  a  Sobolev  embedding 


9 


constant  [33],  and  |MUp(n)  =  (fn(wiwi)p^2)1^p-  We  shall  further  assume  [34]  —  and  verify  a  poste¬ 
riori  —  that  {u(p)  |  p  6  V}  is  a  nonsingular  (isolated)  solution  branch:  (3(u{p))  >  fio  >  0,  V/u  6  V, 
where  fi(z)  =  inf^gx  sup„€X  do(ru,u;2)/||r«||||u||  is  the  inf-sup  parameter  relevant  to  our  nonlinear 
problem.  Numerical  simulations  [29,  32]  demonstrate  that  the  flow  smoothly  evolves  from  a  single¬ 
cell  structure  for  the  lower  Gr  in  V  to  an  inertia-dominated  three-cell  structure  for  the  higher  Gr 
in  V. 

We  may  directly  apply  the  RB  formulation  of  Section  3  to  the  incompressible  Navier-Stokes 
equations  [12,  24,  34].  The  most  significant  new  issue  is  (efficient)  calculation  of  the  nonlinear 
terms.  We  consider  the  inner  Newton  update:  given  a  current  iterate  u//(/u)  =  i^N n(fj) Cn, 
we  must  find  an  increment  6un  G  Wn  such  that  da(SuN,v,u x)  =  V  v  €  W/v;  here 

R(v,p)  =  f{v;p)  -  cl(un(p)}v),  V  v  e  X  is  the  residual.  The  associated  algebraic  equations  are 
thus 

N  N  _ 

{ao(CjiCi)+  5Z  uNn(ll)al{Cjj  Cni Ci)}$uNj 

j~  1  71=  1 

N  N  / 

“  M/o(Ci)  S  2  'U'Nn{p)Q'l{Cj)  Ci)  }UN  Vi  G  N, 

j=l  n-1 

where  we  recall  that  f(v\  p)  =  pf§{v)  and  p  =  Gr. 

We  can  directly  apply  the  offline-online  procedure  described  in  Section  3  for  linear  problems, 
except  now  we  must  perform  summations  both  over  the  affine  parameter  dependence  (rather  trivial 
here)  and  over  the  reduced-basis  coefficients  (of  the  current  Newton  iterate,  ujy(p)).  In  the  online 
stage  —  for  given  new  p  —  at  each  Newton  iteration  u^{p)  — >  Sun  we  first  assemble  the  right-hand 
side  (residual)  —  at  cost  0(iV3);  we  then  form  and  invert  the  left-hand  side  (Jacobian)  —  at  cost 
0(iV3).  The  complexity  of  the  online  stage  is  independent  of  furthermore,  for  our  quadratic 
nonlinearity,  there  is  little  increased  cost  relative  to  the  linear  case.  Unfortunately,  for  a  pth-order 
nonlinearity,  the  online  cost  for  the  residual  assembly  and  Jacobian  formation  will  scale  as  0(iVp+1), 
and  thus  standard  Galerkin  projections  are  viable  only  for  p  =  2  or  at  most  p  —  3  [36].  Fortunately, 
for  larger  p  and  non-polynomial  nonlinearities  —  and  for  non-affine  parameter  dependence  [5]  — 
quite  effective  collocation-like  alternatives  are  available. 

Turning  now  to  a  posteriori  error  estimation,  we  first  “presume”  /?at(^),  a  (to-be-constructed) 
positive  lower  bound  for  the  inf-sup. f parameter  (p)  =  (fO  >  0n(p)  .  >  0,Y/2  k  V. 

We  next  recall  the  dual  norm  ofThe  resfduaf,  £n(p)  =  siipv€^'!R(v;'/ij/llv|] ,  and  introduce  r/v(/i)"  = 
2 p2£N(p)/^(p),  where  p  is  our  L4(Q)-X  embedding  constant.  Finally,  we  define  N*(p)  such  that 
tn(p)  <  1  for  N  >  N*(p);  we  require  N*(p)  <  iVmax,  V/zGD.  (The  latter  is  a  condition  on  Nmax 
that  reflects  both  the  convergence  rate  of  the  RB  approximation  and  the  quality  of  our  inf-sup 
lower  bound.) 

We  may  now  define  our  error  estimator:  for  N  >  N*(p),  A n(p)  =  (j3jj(p)/ p2)(l-  y/l  -  tjv(/z)); 
note  that,  as  £n{p)  — 3 ►  0,  A ^(p)  tends  to  the  “linear  case”  en{p)/Pn{p)‘  Our  main  result  is 
then:  given  any  p  G  V,  for  all  N  >  N*(p),  there  exists  a  unique  (truth  approximation)  solution 
u(p)  G  X  in  the  ball  B(uN{p),/3N{p)/p 2)  =  {z  £  X  \\\z  -  uN(p)\\  <  (3N(p)/p2};  furthermore, 
|| u(p)  -  utf(p) ||  <  Apf(p).  The  proof  [34,  35]  is  a  slight  specialization  of  the  abstract  “Brezzi- 
Rappaz-Raviart”  result  [6,  11];  we  can  further  provide  several  corollaries  related  to  (z)  the  well- 
posedness  of  the  truth  approximation,  and  (ii)  the  effectivity  of  our  error  bound  [34].  (We  may 
also  develop  bounds  for  the  output  of  interest  [34].) 

The  real  challenge,  is  computational:  how  can  we  compute  £n(p)i  P>  an^  Av(aO?  (Note  that, 
armed  with  these  quantities,  we  can  evaluate  rj^(p)  and  hence  verify  N  >  N* (p).)  The  reduced- 
basis  context  is  in  fact  a  rare  opportunity  to  render  the  Brezzi-Rappaz-Raviart  theory  completely 
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quantitative.  To  begin,  we  consider  ew(//):  as  for  the  linear  case,  we  invoke  duality,  our  reduced- 
basis  expansion,  the  affine  parameter  dependence  of  a  (and  /),  and  linear  superposition  to  express 

N  N 

53  Un n([l)  {2/i(C,  £.n)  +  53  UN n'it1)  {2/i(C,  Qnn')  +  {C,n,  £n') 
n= 1  n'—\ 

N  N 

+  X/  Un  n,f  {ft)  Qn'n")  +  X/  (2nn')  2n"nw)j  }  }  >  (4) 

n"= 1  n"'=l 

where  (C,v)  =  /(u),  V  i>  €  AT,  (£n,u)  =  “Oo(Cni«)i  V  v  G  X,  V  n  G  N,  and  (Qnn'>v)  = 

—  iai  (Cn,  Cn'j  v),  V  v  E  X,  V  (n,  n')  G  N2;  the  latter  are  again  simple  (vector)  Poisson  problems. 

We  can  now  readily  adapt  the  offline-online  procedure  developed  in  the  linear  case  [34,  35].  In 
the  online  stage  —  for  each  new  p  —  we  perform  the  sum  (4)  in  terms  of  the  pre- formed  and  stored 
inner  products  (for  example,  (Qnn',  Qn"ri")i  1  <  n,  n",  nm  <  N)  and  the  RB  coefficients  u^n{p), 
1  <  n  <  N  —  at  cost  0(N4).  Although  the  N 4  scaling,  which  arises  due  to  the  trilinear  term  in 
the  residual,  is  certainly  unpleasant,  the  error  bound  is  calculated  only  once:  in  actual  practice,  the 
additional  online  cost  attributable  to  the  dual  norm  of  the  residual  is  not  too  large.  Unfortunately, 
for  a  pth-order  nonlinearity,  the  online  evaluation  of  £w(p)  scales  as  0(N2p),  and  our  approach 
is  thus  viable  only  for  p  =  2.  Fortunately,  for  larger  p  and  non-polynomial  nonlinearities  —  and 
for  non-affine  parameter  dependence  [5]  —  collocation-like  alternatives  are  available;  however,  in 
general,  there  will  be  some  loss  of  rigor . 

We  next  turn  to  the  calculation  of  p.  The  critical  observation  is  that  p  is  the  supremum  of  a 
“Rayleigh-quotient.”  Thus  p  is  related  to  the  smallest  multiplier  of  an  associated  Euler-Lagrange 
nonlinear  eigenproblem  [33]:  (A,  $)  €  (R+,X)  satisfies  —  2A 2 f  V  v  E  X,  for 

ll$lli4(n)  “  the  ground  state  is  denoted  (Amin,^min),  and  p  =  A“*n.  In  practice,  it  may  be 
difficult  to  isolate  the  ground  state;  we  thus  consider  a  homotopy  procedure. 

Towards  that  end,  we  first  introduce  a  parametrized  generalization  of  the  Euler-Lagrange  equa¬ 
tion:  given  a  €  [0,1],  (A(a),^(a))  e  (R+,X)  satisfies  (^,v)  =  2A2(a)[a J + 
(1  —  a) VwGl,  for  normalization  a||'0(a)||£4  +  (1  —  a)||^(a)||22  =  1;  the  ground  state 
is  denoted  (Amin(a),'0min(»)))  and  p  =  A~*n(l).  We  may  now  apply  standard  Newton  continuation 
methods  to  proceed  from  the  known  ground  state  at  a  =  0  —  (Amin(O),<0min(O))  is  thr  lowest 
eiirenpair  of  a  simnle  (vector V  “T^aolnr.i an”  linear  ejgenorobletn  —  to  the  ground  state  bf  interest 
at  a  =  1;  for  sufficiently  small  increments  in  o,  we  will  remain  on  the  desired  (lowest-energy) 
branch.  For  our  particular  domain,  we  find  (offline)  p  =  0.4416;  since  p  is  /i-independent,  no  online 
computation  is  required. 

Finally,  as  regards  the  inf-sup  lower  bound,  /3jv(/z),  we  may  directly  apply  appropriate  exten¬ 
sions  [21,  34]  of  the  procedure  developed  in  Section  5.  The  nonlinear  case  does  present  a  new 
difficulty:  the  parameter  dependence  of  the  (linearized)  operator  is  now  induced  by  the  reduced- 
basis  solution  un(p)  —  in  our  case,  through  the  a\(w ,  u^(/i),  v)  term  * —  and  hence  is  not  known  a 
priori.  Fortunately,  since  un(p)  —>  u(p)  we  may  develop  a  “universal”  lower  bound  for  sufficiently 
large  N ;  the  complications  are  thus  largely  practical  in  nature.  (For  our  particular  problem,  J  =  34 

—  the  sample  is  relatively  small  despite  the  rather  large  range  in  Grashof.) 

We  conclude  with  a  brief  discussion  of  the  adaptive  sampling  procedure  introduced  in  Section 
4.  In  the  nonlinear  case  a  similar  procedure  may  be  pursued,  but  with  two  important  differences. 
First,  as  already  indicated,  (3n{h)  and  hence  /3jv(a0  now  depend  on  the  reduced-basis  solution 
wat(aO;  furthermore,  /?at(/x)  will  only  be  meaningful  for  larger  N.  Thus  in  the  sample  construction 
stage  we  must  replace  Pn{p)  in  A;v(/i)  with  a  simple  but  relevant  surrogate  —  for  example,  a 
piecewise-constant  (over  V)  approximation  to  (3(u(p )).  Second,  in  the  nonlinear  context  our  error 
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Gr  =  l.OEl 

Gr 

=  8.5E4 

N 

TN 

A  N,  rel 

m 

tn 

A  N,  rel 

VN 

3 

2.0  E— 2 

5.0  E-2 

1.0 

OO 

* 

* 

6 

1.2  E-2 

3.0  E-2 

1.0 

2.8  E+l 

* 

* 

9 

4.4  E-3 

1.1  E-2 

1.0 

5.2  E-l 

1.5  E— 4 

14.1 

12 

2.7  E-6 

6.8  E-6 

1.0 

5.8  E-l 

1.7  E-4 

20.5 

15 

3.0  E-7 

7.6  E-7 

1.0 

1.9  E-2 

4.6  E-6 

17.6 

Table  2:  Error  bounds  and  effectivities  for  Gr  =  l.OEl  and  Gr  =  8.5E4. 


bound  is  conditional  —  a  small  solution  to  the  error  equation  is  only  assured  if  Tjy(p)  <  1.  Thus 
the  greedy  procedure  must  first  select  on  argmax^€SF  t/v(m)  —  until  r^(p)  <  1,  V  p  E  EF  * —  and 
only  subsequently  select  on  argmax^^F  A at(/u);  the  resulting  sample  will  ensure  rapid  convergence 
to  a  certifiably  accurate  solution. 

In  Table  2  we  present  AN}Te\(p)  =  and  mil*)  =  &n(p)/ 

\\u(p)  -un(p)\\  818  a  function  of  N  for  p  =  Gr  =  l.OEl  (single-roll)  and  p  =  Gr  =  8.5E4  (three-roll). 
The  indicates  that  N  <  N*(p)  —  tn(p)  >  1:  no  error  bound  is  available.  For  Gr  =  l.OEl,  we 
find  N*(p)  =  1,  and  hence  we  obtain  error  bounds  for  all  AT;  the  error  bound  tends  to  zero  very 
rapidly;  and  the  effectivity  is  0(1)  [21,  34],  For  Gr  =  8.5E4,  we  find  N*(p)  =  9,  and  hence  we 
obtain  error  bounds  only  for  rather  accurate  approximations;  however,  the  error  bound  still  tends 
to  zero  rapidly  with  N  —  our  samples  Sffi  are  constructed  to  provide  uniform  convergence;  and  the 
effectivity  is  still  quite  good.  It  is  perhaps  surprising  that  the  Brezzi-Rappaz-Raviart  theory  —  not 
really  designed  for  quantitative  service  —  indeed  yields  such  sharp  results;  in  fact,  as  ejsr(p)  — *  0, 
the  cruder  bounds  —  in  particular,  p  —  no  longer  play  a  role. 

Finally,  we  note  that  the  online  cost  (on  a  Pentium®  M  1.6GHz  processor)  to  predict  sn{p) 
and  Apj(p)  (and  a  bound  for  the  error  in  the  output,  A^(^)  [34])  is  typically  10ms  and  90ms, 
respectively  —  compared  to  order  minutes  for  direct  finite  element  calculation  of  s(p)  =  £(u(p)). 

7  Parabolic  Equations 

We  consider  here  the  extension  of  the  RB  methods  and  associated  a  posteriori  error  estimators 
described  in  Sections  1-4  to  parabolic  PDEs  —  in  particular,  the  heat  equation;  we  shall  “simply” 
Heat  time  as  an  to  model 

reduction  for  initial- value  problems:  POD  methods  [31];  balanced-truncation  techniques  [20];  and 
even  reduced-basis  approaches  [26].  However,  in  general,  these  frameworks  do  not  accommodate 
parametric  variation  (or,  typically,  rigorous  a  posteriori  error  estimation).)  For  simplicity,  we 
directly  consider  a  X-level  time-discrete  framework  (corresponding  to  Euler  Backward  discretization 
—  we  can  also  readily  treat  higher-order  schemes  such  as  Crank-Nicolson)  associated  to  the  time 
interval  [0 ,£/]:  we  define  T  =  where  tk  =  fcAt,  0  <  k  <  K,  and  At  =  tf/K]  for 

notational  convenience,  we  also  introduce  K=  {1 , . . .  ,if}.  (Clearly,  our  results  must  be  stable  as 
At  — >  0,  K  — >  oo.) 

Given  p  e  V  C  Rp,  we  evaluate  the  (here,  single)  output  s(p,tk)  =  £(u(pytk))y  V k  e  K,  where 
u(p,  tk)  E  X,  V/c  e  K,  satisfies 

A  t_1m(w(^,  tk)  -  u(p,  t^^yv)  +  a(u(py  tk)yV]  p)  =  b(tk)f(v)y  V  v  E  X,  (5) 

with  initial  condition  (say)  u(p,t°)  =  0.  Here  p  and  V  are  the  input  and  input  domain;  u(pytk)y 
V  /c  E  K,  is  our  field  variable;  X  C  Xe  is  our  truth  approximation  subspace  for  Xe  (and  (•,•), 
||  •  ||)  defined  in  Section  2;  a(-,  •;//)  and  m(-,  •)  are  Xe-continuous  and  L2(ft)-continuous  symmetric 
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bilinear  forms,  respectively;  /(•),  £(■)  are  L2(fl)-continuous  linear  forms;  and  b(tk)  is  the  (here, 
single)  “control”  input  at  time  tk. 

We  shall  make  the  following  assumptions.  First,  we  require  that  a  and  m  are  independent  of  time 
—  the  system  is  thus  linear  time-invariant  (LTI).  Second,  we  assume  that  a  and  m  are  coercive: 
0  <  a0  <  a(p)  =  inft)ex[a(v,'y;/x)/||u||2]  and  0  <  a0  <  a(p)  =  irf«eJ*(n)M»htO/IMIia(n)]- 


Third,  we  assume  that  a  depends  affinely  on  p:  a(w,v,p)  =  J2q=i  ag(w,v)  for  q  =  1, . . .  ,Q 

parameter- dependent  functions  ©9(p)  :  V  — »  R  and  parameter- independent  continuous  bilinear 
forms  aq(w,v).  (For  simplicity,  we  assume  that  m,  /,  and  £  are  parameter-independent.) 

To  ensure  rapid  convergence  of  the  reduced-basis  output  approximation  we  shall  need  a  dual  (or 
adjoint)  problem  which  shall  evolve  backward  in  time.  Invoking  the  LTI  property,  we  can  express 
the  adjoint  for  the  output  at  time  tL,  1  <  L  <  K,  as  ipL(p,,tk)  =  \f '(p,tK~L+k),  1  <  fc  <  L\  here 
\k(p,ffc)  G  X,  Vfc  G  K,  satisfies  Af-1m(v,  ^(p,tfc)  -  'S>(fi,tk+1))  +  a(v,  ^(/a,tk);p)  =  0,  Vd  G  X, 
with  final  condition  m(v,  ^(p,  tK+1))  =  £(v),  V  v  G  X. 

We  now  introduce  the  nested  samples  5'^  =  {pPr, . .  ■ ,  P^pt } ,  1  <  Npr  <  -Npr  ,max;  and  S$qu  — 

{p^u, . .  • .  Awdll}>  1  ^  Ndu  <  ^du.max,  where  p  =  (p,ifc)  C  V  =  V  x  T.  Note  the  samples  must  now 
reside  in  the  parameter-time  space  D\  we  also  introduce  separate  (and  different)  samples  for  the 


primal  and  dual  problems.  We  then  define  the  associated  nested  RB  spaces  ^  =  span{£nr  = 
u(pnr  =  (pm  tkn)Pr),  1  <  n  <  Npi},  1  <  Npi  <  NpI>max,  and  —  span{£du  =  ip{fin  — 

(pn,tfcn)du),  1  <  n  <  Ndu},  1  <  Ndu  <  ATdU)max.  Note  that  for  the  primal  basis  we  choose  —  as 
justified  by  the  LTI  hypothesis  —  an  impulse  input,  b(tk)  =  5ik,  V  fc  G  K. 

Our  RB  approximation  is  then:  given  p  G  2>,  evaluate  Sjv(p,  tk)  —  £(un(p,  tk)) 

+  ^/=1i?p^(,^'Ar(p,<^1■_A:+fc,);p,^<::,)  At,  Vfc  G  K,  where  (pr)  uaK/M*)  €  W^r,  VfcsK,  satis¬ 
fies  A t_1m(tijv(p,  tk )  -  tqv(p,  <fc_1),  u)  +  a(uN{p ,  tfc),  u;  p)  =  6(ffc)  /(u), 

Vt)  G  Wjyr,  with  initial  condition  uat(p,  t°)  =  0,  and  (du)  'L/^P,  tk)  G  V  fc  G  K,  sat¬ 

isfies  A t_1m(v,^(p,tfc)  -  #jv(p,<fc+1))  +  o(u,  ^jv(p,tfc);p)  =  0,  V  u  G  W$“u,  with  final  con¬ 
dition  m(u,'F;v(p,<*+1))  =  £(v),  Vug  W$u.  Here,  V  fc  G  K,  flpr(v;p,tfc)  =  b(tk)  f(v)  - 
(At_1m(uAr(p,tfc)  -tijvtp,^-1),^)  +  a(«N(p,<fc),v;p)),  Vu  G  X,  is  the  primal  residual. 

The  offline-online  computational  procedure  is  similar  to  the  elliptic  case  of  Section  3  but  with 
the  added  complexity  of  the  dual  problem  and  the  time  dependence.  In  the  online  stage,  we  first 
APS^rrhle  thejeqffif  ^  RB  “stiffr.cos”  matrices ,-:v  at  cost %lS!SSt^9lySz:^ — 
the  primal  and  dual  problems  —  at  cost  0(N*r  4-  JV|U  4-  K(N%r  4-  Njj);  and  finally  we  evaluate 
the  RB  output  approximation  sN(p]  tk),  Vfc  G  K  —  at  cost  0(KNpiNdu).  The  online  complexity 
is  thus  independent  of  AT,  and  in  fact  not  too  sensitive  (for  our  LTI  system)  to  K. 

We  now  turn  to  a  posteriori  error  estimation.  We  stress  that  the  development  of  the  error 
bounds  is  in  no  way  limited  to  the  RB  approximation  described  here:  we  may  consider  “any” 
stable  ODE  or  PDE  system  and  any  reduced-order  model  To  begin,  we  assume  that  we  are  given 
a(p ;)  :  V  — >  R+  —  a  positive  lower  bound  for  the  coercivity  constant,  a((i)  :  a(p)  >  a(ji)  >  ao  >  0, 

V/i  G  V.  In  our  symmetric  case  a(p)  =  j3(p)  and  thus  a(p)  can  be  constructed  according  to  Section 
4;  in  fact,  thanks  to  coercivity,  much  simpler  procedures  typically  suffice  [27].  We  next  recall  the 
dual  norm  of  the  primal  and  dual  residuals:  V  fc  G  K,  ePNpi{ti,tk)  =  suPt,e*[.RPr(v;p,<fc)/|M|]  and 
e^(p,ffc)  =  sup„ex[i?du(u;p,<fc)/||u||],  where  V  fc  G  K,  i?du(u;p,tfe)  =  -(Ai_1m(u, ^fN{/x,tk) 

—  ^(p,tfc+1))  +  a(v, ^^(p, tfc); p)),  Vv  G  X.  Finally,  we  introduce  the  “spatio-temporal”  energy 
norm,  |||u(p,ffc)|||2  =  m(v(p,tk),v(p,tk))  +  Y%=i  ^  ),  <fc  );  Vv  G  X. 

We  may  now  define  our  error  estimators:  V  fc  G  K,  V  p  G  V,  A^(p,  tk)  =  <5:_2(p)(At^,_1 
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Figure  2:  One  “cell”  of  the  heat  shield. 


SNP>-*fc')2)*;  Awd>,*fc)  =  a-2M(Af  ££=fce$jMfc')2)*;  and 

A S(M*)  =  A £>,*)  Ajfdu(n, tK~k+1). 


(6) 


We  can  then  readily  demonstrate  that  |||u(/z,tfc)— u;v(/i,ffc)|||  <  A^pr(/i,tfc)  and  \s(^,tk)  —  s^(fj,,tk)\ 
<  As(n,  tk),  V  k  eI,V/ieD,V  Wpr  G  Nprimax,  V  Nda  G  Ndu,maj[  —  we  obtain  rigorous  (and,  as  we 
shall  see,  rather  sharp)  upper  bounds  for  the  primal  error,  dual  error,  and  output  error.  (Note  that 
our  particular  form  (6)  assumes  that  \(/(/z,£^+1)  —  here,  /z-independent  —  is  a  member  of 
this  requirement  is  readily  relaxed.) 

The  offline-online  procedure  for  the  computation  of  As(/z,tfc),  V  k  E  K  —  in  particular,  for  the 
calculation  of  the  requisite  primal  and  dual  residual  norms  —  is  similar  to  the  elliptic  case  of  Section 
4  but  with  the  added  complexity  of  the  dual  problem  and  the  time  dependence.  In  particular,  in  the 
online  stage  —  for  any  given  new  \i  —  we  evaluate  the  £^pr  (/z,  tk)2  and  £$*  u  (/z,  tk )2  sums  in  terms  of 
09(/z),  uNn(fi,tk),  ^TVn'QM*)  and  the  precomputed  inner  products  —  at  cost  0(K(N2T+N$u)Q2). 
Thus,  all  online  calculations  are  indeed  independent  of  AT. 

We  now  turn  to  a  particular  numerical  example.  We  consider  the  design  of  a  heat  shield  (one 
cell  of  which  is  shown  in  Figure  2):  the  left  boundary  ddout  is  exposed  to  a  temperature  unity 
and  Biot  number  Bi0ut  “source”  for  t  E  [0,  £/];  the  right  boundary  as  well  as  the  top  and  bottom 
(symmetry)  boundaries  are  insulated;  and  the  internal  boundaries  dflm  —  corresponding  to  three 
square  cooling  ch&r<~  Js  Hie .exposed::to:''a:t'emp.erature  zero  ^and  B;ot  number  Biin  ;“sink.”  ,:Our 
input  parameter  is  hence  /z  =  (m(i)iM(2))  —  (Bi0ut5Biin)  =  [0.01,0.5]  x  [0.001,0.1];  our  output 
is  the  average  temperature  of  the  structure  —  a  surrogate  for  the  maximum  temperature  of  the 
(to-be-protected)  right  boundary  for  t  E  [0,  oo[ . 

The  underlying  PDE  is  the  heat  equation.  The  (appropriately  non-dimensionalized)  govern¬ 
ing  equation  for  the  temperature  u(p^tk)  E  X  is  thus  (5),  where  X  is  a  linear  finite  element 
truth  approximation  subspace  (of  dimension  (exploiting  symmetry)  A f  =  1,396)  of  Xe  =  J?1(fi); 
o(u>,  V\n)  =  fn  VwVv+w  fnoutwv+m  fninwv>  m(w,v)  =  Jnwv;  f(v,n)  =  /i(i)  fQout  v,  which 
is  now  (affinely)  parameter-dependent;  b(tk)  =  1,  V/c  G  K;  and  (w,  v)  =  fn  Vw-Vv+0.01  Janout  wv+ 
0.001  Jqq,  wv  —  hence  we  may  choose  a(fi)  =  1.  The  output  is  given  by  s(fi,  tk)  —  t(u(fi,tk)), 
where  i(v)  =  |f2|-1  fnv. 

We  now  present  numerical  results.  Our  “optimal”  primal  and  dual  samples  are  constructed 
(separately)  by  procedures  similar  to  the  greedy  approach  described  for  the  elliptic  case  in  Section 
4:  at  each  step  (say,  for  the  primal)  we  select  the  parameter  value  n*  for  which  A ^prGMK)  is 
maximized;  we  then  select  the  time  tk *  for  which  e^r  r(^Vfc)  is  maximized.  In  Table  3  we  present, 
as  a  function  of  Npi  (=  ATdu),  A^  ^,  ifjpr,  AJ^/and  fjs:  Aj^,  is  the  maximum  over  STest  of 
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Npr 

Apr  , 

max,  r  el 

fjpi 

As  , 

max,rel 

—5 

V 

4 

1.6E-00 

5.44 

1.6E-00 

95.63 

8 

6.3E-02 

1.55 

6.7E-03 

30.92 

12 

1.0E-02 

1.03 

2.6E-04 

8.43 

16 

3.2E-03 

1.02 

1.5E-05 

11.45 

20 

8.8E-04 

1.01 

1.1E-06 

17.43 

Table  3:  Convergence  results  for  the  heat  equation. 

&Hpr{v,tK)/\\\uN{fiu,tK)\\\,  rjpt  is  the  average  over  HTest  xT  of  Ap^pT{fi,tk)/\\\u(fj,,tk)-uN(fJ,,tk)\\\, 
A^lax  r(:i  is  the  maximum  over  Hxest  of  Asn(/j, tK) /\s^(f.is,  tK)\,  and  rjs  is  the  average  over  STest  of 
AaN(n,  tv{fj,))/\s{n,  m/j,))  -sn{ij.,  tn{fi))\.  Here  “Test  €  (£>)400  is  a  random  input  sample  of  size  400; 
Hu  =  argmax^€HTest  |||u;vinax(/i,f*')|||,  A*.  -  arg  max^gHTest  |«]w(a*.**)I  (note  the  output  grows 
with  time),  and  tr)(M)  =  argmaxt*6T  |s(^,ffc)  —  s^(^t,tfc)|. 

Finally,  we  note  that  the  calculation  of  tk)  and  A sN(/j,,tk),  V  k  €  K,  is  (say,  for  iVpr  = 
JVdu  =  12)  roughly  40  x  faster  than  direct  calculation  of  the  truth  approximation  output  s(/x,ffc )  = 
e(u(ii,tk)),  Vfc  €  K.  We  may  thus  work  with  sN(fi,  tk)  +  AsN(fi,tk)  as  a  certifiably  conservative 
(upper  bound)  and  accurate  surrogate  for  the  average  temperature  s(fi,tk)  in  truly  interactive 
design  exercises. 
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Project  3  —  Computing  Bounds  for  Functional  Output  of  Solutions  of  PDE’s 
Project  3(a): 


Computing  Bounds  for  Functional  Output  of  Solutions  of  PDE’s 


Investigator:  Jaime  Peraire 

Department  of  Aeronautics  and  Astronautics 
Room  37-451 

Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139 

Background:  Exact  Bounds  and  Certificates 

Existing  techniques  for  approximating  the  solutions  of  PDE’s  rely  on  the  experience  of  the  user 
to  estimate  a-priori  the  size  of  the  discretization  required  to  resolve  the  various  problem  features. 
Failure  to  appropriately  do  so,  may  result  in  results  which  are  either  too  expensive  to  obtain  or, 
worse  yet,  uncertain.  Modern  a-posteriori  and  mesh  adaptivity  methods  have  the  potential  to 
alleviate  this  problem,  but  not  to  eliminate  it  completely,  i.e.  a  saturation  hypothesis  needs  to  be 
made  that  can  not  be  verified  a-priori. 

For  a  number  of  years  we  have  been  engaged  in  the  development  of  techniques  for  computing 
strict  bounds  for  functional  outputs  of  the  exact  solution  of  partial  differential  equations.  Our 
approach  draws  on  previous  work  on  a-posterior  error  estimation  and  complementary  energy  ideas 
well  known  in  the  mechanics  community.  We  have  generalized  the  previous  work  in  a  number  of 
ways.  Relative  to  existing  a-posteriori  error  estimation  techniques,  our  approach  delivers  bounds 
with  respect  to  the  exact  solution  of  the  PDE  as  opposed  to  just  an  estimate  of  the  computed  error. 
Relative  to  the  complementary  energy  methods,  which  are  only  capable  of  delivering  bounds  for 
the  energy  in  problems  for  which  a  variational  principle  exists,  we  are  able  to  handle  general  linear 
outputs,  some  non-linear  outputs,  and  some  non-symmetric  problems  both  linear  and  non-linear 
for  which  variational  principles  may  not  exist.  ^ 

The  startipp*  point  for  bounds  .prqcedure-is^.^^^elfihenf.;app^qjiip;fti5QP ^ojthe  sohition 

and  to  the  output  dependent  adjoint  solution.  These  approximations  are  then  post-processed  to 
yield  the  so  called  inter-element  hybrid  fluxes.  The  hybrid  fluxes  are  then  used  as  data  for  the 
computation  of  locally  equilibrated  stress  fields.  The  final  expression  for  the  bounds  is  obtained 
by  calculating  appropriate  norms  of  the  stress  fields.  It  is  shown  that  the  computed  bounds  are 
uniformly  valid  regardless  of  the  size  of  the  underlying  coarse  discretization,  but  as  expected,  their 
sharpness  depends  on  the  accuracy  of  the  approximated  solutions.  A  mesh  adaptive  procedure  has 
also  been  developed  which  can  be  used  to  determine  the  bounds  to  a  preset  level  of  accuracy. 

An  attractive  feature  of  our  approach  is  that  the  piecewise  polynomial  equilibrated  stress-like 
fields,  which  are  computed  as  part  of  the  bound  process,  can  be  used  as  certificates  to  guarantee 
the  correctness  of  the  computed  bounds.  It  turns  out  that  given  a  stress  field,  it  is  easy  to  check 
whether  this  field  corresponds  to  a  valid  certificate,  and  in  the  affirmative  case,  it  is  straightforward 
to  determine  the  value  of  the  output  that  it  can  certify.  In  particular,  the  stress  fields  need  to  satisfy 
continuity  of  normal  tractions  across  elements,  and  membership  of  an  appropriate  space. 

The  idea  of  a  certificate  that  is  computed  simultaneously  with  the  solution  has  many  attractive 
features.  In  particular,  a  certificate  consisting  of  the  data  set  necessary  to  describe  the  piecewise 
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polynomial  stress-like  fields  could  be  used  to  document  the  computed  results.  We  note  that  exer¬ 
cising  the  certificate  does  not  require  access  to  the  code  used  to  compute  it,  and  can  be  done  with 
a  simple  algorithm  which  does  not  require  solving  a  system  of  equations.  A  very  important  point 
is  that,  if  a  certificate  meets  all  the  necessary  conditions,  which  in  turn  are  easy  to  verify,  then, 
there  is  no  need  to  certify  the  code  used  to  compute  it.  In  practice,  the  size  of  these  certificates 
depends  on  the  required  level  of  certainty.  As  expected,  we  shall  find  that  high  levels  of  certainty, 
i.e.  small  bound  gaps,  will  often  require  longer  certificates  (larger  data  sets)  than  those  required 
to  certify  less  sharp  claims. 

To  date,  we  have  been  successful  at  developing  algorithms  for  computing  bounds  for: 

-  Linear  Functional  Outputs  of  the  Convection-Diffusion-Reaction  Equation 

-  Linear  Functional  Outputs  of  the  Linear  Elasticity  Equations 

-  Energy  Release  Rates  (J-Integral)  in  Linear  Elasticity 

-  Linear  Functional  Outputs  for  the  Stokes  Equation 


1  Particular  Achievements 

The  focus  of  this  project  has  been  the  extension  of  our  methodology  to  the  computation  of 
strict  upper  and  lower  bounds  for  for  the  collapse  load  in  limit  analysis. 

Limit  analysis  is  relevant  in  many  practical  engineering  areas  such  as  the  design  of  mechanical 
structures  or  the  analysis  of  pressurized  vessels.  Whereas  linear  elastic  analyses  are  typically  used 
to  determine  the  performance  of  the  structure  (usually  characterized  by  deflections)  under  the 
so  called  service  loads,  limit  analysis  is  used  to  determine  the  collapse  load.  Assuming  a  rigid, 
perfectly-plastic  solid  subject  to  a  static  load  distribution,  the  problem  of  limit  analysis  consists  of 
finding  the  minimum  multiple  of  this  load  distribution  that  will  cause  the  body  to  collapse.  This 
collapse  multiplier  results  from  solving  an  infinite  dimensional  saddle  point  problem,  where  the 
internal  work  rate  is  maximized  over  an  admissible  set  of  stresses  -defined  by  a  yield  condition-  and 
minimized  over  the  linear  space  of  kinematically  admissible  velocities  for  which  the  external  work 
rate  equals  the  unity.  When  strong  duality  is  applied  to  this  saddle  point  problem,  the  well-known 
crrver.  (and  tequivalerit)_atatic  and.kiaematic^rine1pler5.x>f  limit  analysis  arise  in  this^prpject,  ap  • 
efficient  procedure  to  compute  strict  upper  and  lower  bounds  for  the  exact  collapse  multiplier  has 
been  developed,  with  a  formulation  that  explicitly  considers  the  exact  convex  yield  condition.  The 
approach  consists  of  two  main  steps.  First,  the  continuous  problem,  under  the  form  of  the  static 
principle,  is  discretized  twice  (one  per  bound)  by  means  of  different  combinations  of  finite  element 
spaces  for  the  stresses  and  velocities.  For  each  discretization,  the  interpolation  spaces  are  chosen  so 
that  the  attainment  of  an  upper  or  a  lower  bound  is  guaranteed.  The  second  step  consists  of  solving 
the  resulting  discrete  nonlinear  optimization  problems.  Towards  this  end,  they  are  reformulated  into 
the  canonical  form  of  Second-order  Cone  Programs,  which  allows  for  the  use  of  primal-dual  interior 
point  methods  that  optimally  exploit  the  convexity  and  duality  properties  of  the  limit  analysis 
model  and  guarantee  global  convergence  to  the  optimal  solutions.  To  exploit  the  fact  that  collapse 
mechanisms  are  typically  highly  localized,  a  novel  method  for  adaptive  meshing  is  introduced  based 
on  local  bound  gap  measures  and  not  on  heuristic  estimates.  The  method  decomposes  the  total 
bound  gap  as  the  sum  of  positive  elemental  contributions  from  each  element  in  the  mesh,  and 
refines  only  those  elements  which  are  responsible  for  the  majority  of  the  numerical  error.  Finally, 
stand-alone  computational  certificates  that  allow  the  bounds  to  be  verified  independently,  without 
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recourse  to  the  original  computer  program,  are  also  provided.  This  removes  the  uncertainty  about 
the  reliability  of  the  results,  which  frequently  undermines  the  utility  of  computational  simulations. 

1.1  Application  Example 

For  illustration  purposes,  we  show  two  computations  on  the  same  geometry.  One  is  carried  out 
using  a  linear  elastic  model  and  the  other  one  carried  out  using  limit  state  analysis.  The  figure 
below  shows  a  square  thin  plate  with  two  rectangular  cut-outs.  The  plate  is  loaded  with  a  uniform 
traction  on  its  boundary,  and  as  a  consequence,  the  material  develops  a  strain  induced  stresses. 
Due  to  symmetry  only  one  quarter  of  the  plate  is  analyzed. 


r° 


r3 


Figure  1:  Thin  square  plate  with  rectangular  holes  under  a  plane  stress  loading. 

1.1.1  Linear  Elastic  Analysis 

In  this  case,  we  are  interested  in  the  maximum  displacement  that  the  portion  of  the  boundary 
T°  would  experience  as  a  result  of  the  loading.  Therefore,  our  output  of  interest  becomes  the 
average  horizontal  displacement  over  T®.  The  figure  below  shows  the  sequence  of  adapted  meshes 
automatically  generated  during  this  computation.  The  table  shows  the  lower  and  upper  bounds 
computed,  s±,  as  well  as  the  bound  gap,  A,  as  a  function  of  the  number  of  elements  in  the  mesh. 
We  point  out  that  this  bounds  are  strict  with  respect  to  the  analytical  solution  of  the  problem  and 
that  a  certificate  for  these computations  could  fee  ehSily  gelierated.  —  - 
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1966 

2532 

3069 

3564 

A 

.1074 

.1821 

.1217 

.0719 

.0375 

.0242 

.0157 

.0117 

.0083 

s+ 

S- 

.3695 

.2620 

.3894 

.2072 

.3688 

.2470 

.3508 

.2789 

.3375 

.2999 

.3339 

.3096 

.3292 

.3134 

.3282 

.3165 

.3262 

.3179 

Table  1:  Computed  bounds  in  a  series  of  adaptively  /i-refined  meshes. 


1.1.2  Limit  Analysis 

Here,  the  problem  is  studied  using  the  limit  analysis  theory.  Given  a  yield  criterion  for  the 
material  (here  the  von  Mises  yield  criterion  has  been  used),  we  want  to  calculate  the  maximum 
load  multiplier  that  can  be  applied  to  the  plate  before  the  plate  collapses.  The  problem  has  been 
solved  using  the  developed  formulation  and  the  SDP3T  second  order  cone  solver.  It  is  worth  noting 
that  in  this  case  the  collapse  mechanism  consists  of  two  localized  shear  bands  which  concentrate 
most  of  the  deformation  whereas  the  rest  of  the  plate  behaves  very  much  like  a  rigid  body.  Below  we 
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Figure  2:  Sequence  of  adapted  meshes  for  the  average  displacement  output  over  T°.  Number  of 
elements:  108,165,280,405  and  538 


illustrate  the  computed  results.  The  left  figure  shows  the  contours  of  equal  horizontal  displacement, 
the  middle  figure  shows  the  contours  of  equal  effective  deformation,  and  the  right  figure  shows  the 
convergence  of  the  upper  and  lower  bounds  for  the  load  multiplier,  as  the  grid  is  refined. 


Figure  3:  Contours  of  equal  horizontal  displacement,  contours  of  equal  effective  deformation,  and 
convergence  of  the  upper  and  lower  bounds  for  the  load  multiplier  as  the  grid  is  refined. 


2  Personnel  Supported 

The  funds  provided  for  this  project  were  used  to  partially  support  H.  Ciria  during  his  S.M. 
degree.  The  title  of  his  thesis  is  “Computation  of  Upper  and  Lower  Bounds  in  Limit  Analysis 
Using  Second  Order  Cone  Programming  and  Mesh  Adaptivity” 

3  Interactions 

The  work  reported  here  involves  numerous  collaborations  with  other  research  groups  as  well  as 
institutions.  The  idea  of  certification  of  computational  results  is  of  direct  interest  to  the  research 
group  at  Sandia  National  Laboratories  with  which  we  have  a  longstanding  collaboration.  The  work 
developed  under  this  project  has  been  presented  at  conferences  and  meetings  over  the  last  year 
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and  the  initial  feedback  has  been  positive.  At  the  personal  level,  we  have  interacted  closely  with 
Professors  A.T.  Patera  and  P.  Parrilo  from  MIT,  A.  Huerta  and  N.  Pares  from  UPC,  Barcelona,  J. 
Bonet  from  UC  Swansea,  and  K.C.  Toh  from  NUS. 

4  Publications  Resulting  from  this  Research 

The  following  papers  have  resulted  directly  from  this  research: 

•  H.  Ciria  and  J.  Peraire,  Computation  of  Upper  and  Lower  Bounds  in  Limit  State  Analysis 
using  Second-Order  Cone  Programming  and  Mesh  Adpativity,  presented  at  the  9th  ASCE 
Speciality  Conference  on  Probabilistic  Mechanics  and  Structural  Reliability,  June  2004,  Al¬ 
buquerque. 

•  H.  Ciria  and  J.  Peraire,  Upper  and  Lower  Bounds  in  Limit  State  Analysis  using  a  Second- 
Order  Cone  Programming  Formulation,  in  preparation. 

•  H.  Ciria,  J.  Bonet  and  J.  Peraire,  Mesh  Adaptive  Algorithm  for  Limit  State  Analysis,  in 
preparation 

The  following  papers  are  related  to  this  research.  They  present  rigorous  upper  and  lower  bound 
'algorithms  for  different  classes  of  equations: 

•  Sauer-Budge  A.M.,  Bonet  J.,  Huerta  A.,  Peraire  J.,  Computing  bounds  for  linear  functionals 
of  exact  weak  solutions  to  Poisson’s  equation,  SIAM  Journal  on  Numerical  Analysis ,  Vol  42, 
4,  1610-1630,  2004. 

•  Sauer-Budge  A.M.,  Peraire  J.,  Computing  bounds  for  linear  functionals  of  exact  weak  solu¬ 
tions  to  the  advection-diffusion-reaction  equation,  SIAM  Journal  on  Scientific  Computing , 
Vol  26,  2,  636-652,  2004. 

•  Pares  N.,  Bonet  J.,  Huerta  A.,  Peraire  J.  The  computation  of  bounds  for  linear-functional 
outputs  of  weak  solutions  to  the  two-dimensional  elasticity  equations,  submitted  to  Comp. 
Meths.  Appl.  Mech.  Engr.  2004. 

c  Xaan  Z.C,  Lee  RvH.,  Patera  for-dhe* 

J-integral  in  two-dimensional  linear  elasticity,  SMA  Symposium,  January  2004. 

•  Xuan  Z.C.,  Pares  N.,  Peraire  J.,  Computing  upper  and  lower  bounds  for  energy  release  rates 
in  linear  elasticity,  submitted  to  Comp.  Meths.  Appl.  Mech.  Engr.,  2004. 


5  Future  Work 

We  feel  that  this  research  can  be  extended  in  a  number  of  meaningful  ways.  In  particular,  we  are 
looking  into  deformation-plasticity  constitutive  models  which  can  be  interpreted  as  a  generalization 
of  the  linear  elastic  and  limit  analysis  theories.  To  our  knowledge,  no  methods  have  been  proposed 
which  are  able  to  deliver  bounds  for  these  types  of  theories. 
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Project  3(b): 


Bounds  on  Functional  Outputs  and  Semidefinite  Programming2 


Investigator:  John  C.  Doyle 

Control  and  Dynamical  System  107-81 
1200  E.  California  Blvd. 

California  Institute  of  Technology 
Pasadena,  CA  91125 

Our  work  in  this  area  has  concentrated  on  the  algorithmic  construction  of  Lyapunov-like  func¬ 
tions  that  certify  upper  and  lower  bounds  on  functional  outputs  of  PDEs.  These  functions  are  a 
generalization  of  barrier  functions  for  ODEs  [5],  a  framework  that  can  be  used  to  answer  many 
analysis-type  questions  in  systems  theory,  such  as  model  invalidation  based  on  data  and  hybrid 
system  verification  [6].  Apart  from  the  problem  of  estimating  functional  outputs,  we  have  also 
tackled  the  problem  of  assessing  stability  of  steady-states  of  parabolic  PDEs  through  the  construc¬ 
tion  of  Lyapunov  functions.  Both  results  have  shown  that  semidefinite  programming  in  general  and 
sum  of  squares  in  particular  open  new,  unexplored  possibilities  and  provide  a  unique  perspective 
for  future  research  in  systems  described  by  Partial  Differential  Equations.  A  technical  paper  is 
currently  under  preparation. 

1  Estimating  bounds  for  functional  outputs 

We  have  developed  techniques  for  computing  upper  and  lower  bounds  on  functional  outputs  for 
all  three  types  of  PDEs  —  elliptic,  parabolic  and  hyperbolic  —  in  their  standard  form  with  forcing 
terms  and  the  necessary  boundary  and  initial  conditions.  The  functional  outputs  may  be  spatial 
and/or  temporal  averages  of  the  state  or  functions  of  state  of  the  system. 

In  particular  we  developed  two  methodologies  that  are  complementing  each  other  in  many 
respects.  First,  for  elliptic  and  parabolic  type  PDEs  we  use  the  maximum  principle  to  estimate 
functional  outputs.  For  PDEs  where  the  maximum  principle  does  not  hold,  barrier  certificates  can 
be  used  instead.  Here  we  outline  these  two  techniques  through  illustrative  examples.  In  a  later 
section,  we  concenferate::bn'uthe..issiie»-6£..assessing  stability  of  steadyrStatea;foriPDEs.,  througfr.theJ 
construction  of  Lyapunov-type  certificates. 


1.1  PDEs  with  a  ‘maximum’  principle 

This  method  resembles  the  method  of  subharmonics  in  Perron’s  process  [3],  and  is  valid  for 
elliptic  and  parabolic  PDEs.  It  is  based  on  the  construction  of  approximations  to  the  true  solution. 
To  illustrate  the  technique,  consider  a  1-D  Poisson  equation  defined  in  a  region  fi  =  [0, 1]: 


d2u  _ 
dt2  ~ 


u(0)  =  A0, 


u(l)  =  A\. 


Suppose  the  output  functional  y  =  J’q  udt  needs  to  be  estimated.  To  do  this,  we  search  for  a 
polynomial  function  B(t)  so  that 


d?B 

dt2 


-f<  o, 


-7  <  0, 


B(  0)  =  Aq, 


B(  1)  =  Ax 


2Written  by  Antonis  Papachristodoulou.  Email:  antonisOcds.caltech.edu 
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The  first  condition  guarantees  that  B(t)  >  u(t );  the  second  condition  imposes  f*  u(t)dt  <  B(t)dt  < 
7  and  therefore  a  good  upper  bound  can  be  obtained  as  the  order  of  B  is  increased  and  7  is  op¬ 
timized.  Lower  bounds  can  also  be  obtained  by  changing  the  direction  of  the  inequalities  in  the 
above  conditions. 

The  non-negativity  conditions  that  the  function  B  has  to  satisfy  can  be  relaxed  to  the  existence 
of  a  sum  of  squares  decomposition;  using  SOSTOOLS  [7],  such  functions  can  be  constructed  and 
upper  and  lower  bounds  on  y  can  be  estimated.  The  same  technique  can  be  used  for  parabolic 
equations  which  satisfy  the  maximum  principle. 

The  disadvantage  of  this  technique  is  that  the  output  functional  has  to  be  linear  in  the  unknown 
coefficients  of  u ,  for  otherwise  the  problem  is  not  convex  anymore  and  semidefinite  programming 
cannot  be  used  directly  to  construct  the  functions  B.  Moreover  for  hyperbolic  equations,  this  tech¬ 
nique  cannot  be  used.  To  tackle  these  two  cases,  we  introduced  the  method  of  barrier  certificates. 


1.2  Using  barriers  to  the  solution 

This  methodology  centers  on  disproving  that  certain  bounds  can  be  obtained,  by  constructing 
appropriate  state-barrier  functions.  More  theory  can  be  found  in  [5],  where  the  problem  of  model 
invalidation  based  on  data  is  tackled.  In  particular,  consider  a  system  of  the  form: 

^  =  f(x,t),  x(to)  =  x0  (1) 

where  x  €  X  C  Rn  and  /  is  such  that  solutions  exist  and  are  unique  in  X.  Suppose  we  wanted 
to  disprove  that  at  time  tf,  the  solution  Xf  €  Xf  can  never  be  achieved.  This  can  be  done  by 
exhibiting  a  function  B(x,t )  with  the  following  properties: 


B(xf,tf)  -  B(xo,to)  >  e, 


dB 

dt 


dB(x,t )  dB(x,t ) 


+ 


f(x,t )  <  0, 


for  xj  £  Xf}  e  >  0 
for  t  £  [to,  ^i],  x  £  X. 


dt  dx 

The  first  condition  ensures  that  B  at  the  final  time  has  a  value  that  is  strictly  greater  than  at  the 
initial  time,  whereas  the  second  condition  ensures  the  opposite:  that  B  decreases  as  the  system 
evolves.  Therefore  the  set  Xj  can  never  be  reached.  To  use  the  above  result  algorithmically,  we 
first  capture  the  sets  X  and  Xf  by  a  set  of  inequalities  To  verify  the  above  conditions,  we  use 
two  tools:  the  S-proeedure  (a  generalization  of  the  Lagrange  multiplier  method  in  optimization) 
to  adjoin  these  inequalities  to  the  corfespondiiig  conditions;  aridThe'surii  of  squares  decomposition 
and  SOSTOOLS  to  construct  B  algorithmically. 

The  above  methodology  has  been  fully  developed  for  systems  described  by  Ordinary  Differential 
Equations.  In  this  section  we  will  use  this  technique  to  estimate  functional  outputs  with  more 
general  kernels,  and  treat  the  wave  equation.  Consider  first  the  one-dimensional  Poisson’s  equation 

d2u  1  u(t0)  =  Ao,  u(t  0  =  Ax. 


dt2 


f  for  t  £  [toxhi 


Let  the  functional  output  to  be  estimated  be 

rh 


m  =  L\  7t 


2  fu  I  dt 


It  is  easy  to  verify  that  the  above  setup  is  equivalent  to 


d_ 

dt 


“  Ui  “ 

U2 

U2 

= 

f 

.  U3  . 

_  -  2/ui  _ 

,  y  =  u3(ti),  ui(t0)  =  Ao,  u3(t0)  =  0,  ux(ti)  =  Ai. 
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We  have  thus  turned  our  problem  into  the  aforementioned  standard  form.  Now,  we  can  postulate 
appropriate  regions  Xf  to  be  disproved,  yielding  upper  and  lower  bounds  on  ^3(^1). 

The  case  of  the  1-D  wave  equation  can  be  treated  in  the  same  way.  Let 

^  +  k2u  =  /,  t  €  [to,  ti],  '  u(fo)  =  Ao,  u{t  1)  =  A\ 

Let  the  output  condition  be 


The  above  is  equivalent  to 


"  Ui  “ 

U2 

U2 

= 

-k2ui  +  / 

m  W3  _ 

Ul 

y  =  us(ti),  ui{tQ)-ao,  us(t0)  =  0,  ui(ti)  =  ai. 


Again  the  system  is  in  the  standard  form  and  upper  and  lower  bounds  can  be  estimated  directly. 

The  question  that  emerges  is  whether  PDEs  in  more  variables  can  be  treated.  To  do  this, 
we  need  to  generalize  the  notion  of  barrier  functions  to  the  notion  of  barrier  functionals.  The 
complication  is  superficial;  for  example,  consider  the  one-way  wave  equation 


ut  =  ux  +  tx(7t  —  x),  ii(0,£)  =  0,  u(tTi  t)  =  0,  u{x,  0)  —  0. 


Suppose  we  wanted  to  estimate: 


u(0.5)  =  [  u(x,0.5)dx 

Jo 

Here,  we  can  construct  a  functional  of  the  form: 


B(v,  t)  =  ao(t)  +  ai(t)  [  udx  +  a2(t) 
Jo 


(2) 


Now  we  have: 

dB  dB  r7r  4  f 77  r 

—  =  —  +  ai(t)  I  utdx  +  2a%(t)  /  udx  /  utdx  +  ... 

Prom  the  describing  equation,  through  integration  by  parts  we  get  utdx  =  The  barrier 
methodology  can  be  applied  to  this  case  directly,  through  which  we  can  get  upper  and  lower 
bounds  on  u(x,0.5)dx  that  converge  to  the  true  value  of  0.6296. 

The  same  technique  can  be  used  to  analyze  other  types  of  PDEs,  even  nonlinear.  Consider  for 
example  Burger’s  equation: 


Uf  -f-  1 Lltx  —  w(0,  0  —  ^(1,  t )  —  0,  0)  1  X. 

Suppose  we  were  interested  in  a  bound  on  the  value  of  at  x  —  1  at  a  specified  or  at  an  infinite 
time  horizon.  Consider  a  function  of  the  form  (2).  Then  it  is  easy  to  see  that  now 


utdx  =  vux  |J  — 


Using  B  as  a  barrier  function,  we  can  estimate  the  value  of  vux  |J.  With  our  data  this  value  is  0.5, 
which  is  what  the  upper  and  lower  values  obtained  using  the  barrier  method  give. 
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2  Stability  analysis  of  systems  described  by  PDEs 

Our  work  on  the  algorithmic  construction  of  barrier  functionals  to  estimate  functional  outputs 
for  PDEs  was  motivated  by  the  related  issue  of  stability  analysis  of  infinite  dimensional  systems  [4, 
2,  1].  For  this,  we  used  a  methodology  based  on  the  construction  of  functional  certificates  of  integral 
type  whose  kernels  are  polynomials.  The  construction  is  done  using  semidefinite  programming  and 
SOSTOOLS. 

Stability  notions  for  systems  described  by  PDEs  can  be  found  in  [2] .  Through  a  careful  definition 
of  a  nonlinear  system  (semigroup)  S(t)  on  a  complete  metric  space  C  with  metric  p  and  the  orbit 
through  a  point  x  in  this  space,  one  defines  what  is  an  equilibrium  point  Stability  is  then  defined 
through  an  e-8  argument,  under  the  above  metric,  and  to  ensure  it,  one  can  construct  a  Lyapunov 
function,  i.e.  continuous  real-valued  function  V  on  C  such  that 

V(x)  4  {V(S(t)x)  -  V(x)}  <  0 

for  all  x  eC,  not  excluding  the  possibility  V(x)  =  -oo.  The  Lyapunov  stability  theorem  is  similar 
to  the  one  in  finite  dimensions,  and  reads 

Theorem  1  (Lyapunov  theorem)  [2]  Let  {S(t),t  >  0}  be  a  dynamical  system  on  C}  and  let  0  be  an 
equilibrium  point  in  C.  Suppose  V  is  a  Lyapunov  function  which  satisfies  V(0)  —  0,  V{x)  >  c(||a;||) 
for  x  G  C,  ||rr ||  =  p{a;,0},  where  c(-)  is  a  continuous  strictly  increasing  function ,  c(0)  =  0  and 
c(r )  >  0  for  r  >  0.  Then  0  is  stable.  If  in  addition  V(x)  <  — ci(||rc||)  where  ci(-)  is  also  continuous , 
increasing  and  positive  with  c\  (0)  =  0,  then  0  is  asymptotically  stable. 

Although  there  is  a  converse  theorem,  finding  a  Lyapunov  function  in  any  particular  case  is  usually 
difficult,  especially  finding  one  that  satisfies  all  the  above  requirements.  Instead,  the  sum  of  squares 
decomposition  can  be  used  to  construct  them  algorithmically. 

Consider  for  example  the  system 

ut  =  uxx  +  pu3  +  au  (3) 

tz(±l,i)  =  0 

with  p  <  0.  From  the  linearization  of  this  system  about  the  zero  steady  state,  stability  is  guaranteed 
Zoi  a  <  ^i?/4.  Ilere  we- are4nteresteck4i^on^tructing*4ie^3yapunov-fu'net:orji  - 

V(u)  =  j  a(u{rj,t),r])dri  (4) 


where  a  is  a  polynomial  in  its  arguments,  of  order  at  least  2  in  u.  When  V  is  considered,  the 
system  dynamics  come  into  play  and  through  integration  by  parts  of  the  resulting  expression,  the 
boundary  conditions  appear.  In  particular,  the  following  ‘by-product’  conditions  can  be  obtained 


+  u\ 


>d2a 

W 


u 


d3a 

_du2drjUrl 


d3a 

duvd7]2 


} 


dp  =  0 


as  a  result  of  integration  by  parts  of  the  term  ^u^drj.  Also,  we  have  used  the  fact  that 
=  0  and  the  fact  that  =  0  since  a  is  a  polynomial  of  degree  at  least  2  in  u. 

Notice  that  the  following  is  also  true: 


J  r]nu2dr}  =  -  J  nr)nu2dr)  -2  J  rf^uu^dr] 


27 


i.e. 


J  [(1  +  ri)r)nu2  +  2?7n+1tmr)]  drj  =  0 


Stability  in  L2  can  then  be  guaranteed  by  constructing  a  polynomial  a(u,  rj)  and  a  positive  definite 
function  <p  such  that  the  following  hold: 

1.  a(u,  rj)  -  <p(u)  >  0  when  7?  €  [-1, 1], 

2.  -a^u-u2^+u  un  +  ggfa] - ew 2  >  0  when  V  €  [-1, 1]  and  (l+n^u^^u^  = 
0  is  satisfied. 

This  is  because  the  two  conditions  guarantee  the  existence  of  a  Lyapunov  functional,  given  by  (4), 
that  is  positive  definite  and  whose  derivative  is  negative  definite.  Numerical  tests  show  that  as  the 
order  of  a  with  respect  to  rj  is  increased,  the  conservativeness  in  the  bounds  obtained  is  reduced. 
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